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ABSTRACT 

Starting from the hypothesis that the Galaxy’s dark halo responded adiabatically to 
the infall of baryons, we have constructed a self-consistent dynamical model of the 
Galaxy that satishes a large number of observations, including measurements of gas 
terminal velocities and masers, the kinematics of a 180 000 giant stars from the RAVE 
survey, and star count data from the SDSS. The stellar disc and the dark halo are both 
specified by distribution functions (dfs) of the action integrals. The model is obtained 
by extending the work of Piffl, Penoyre & Binney from the construction of a single 
model to a systematic search of model space. Whereas the model of Piffl, Penoyre & 
Binney violated constraints on the terminal-velocity curve, our model respects these 
constraints by adopting a long scale length = 3.66 kpc for the thin and thick discs. 
The model is, however, inconsistent with the measured optical depth for microlensing 
of bulge stars because it attributes too large a fraction of the density at i? < 3 kpc 
to dark matter rather than stars. Moreover, it now seems likely that the thick disc’s 
scale-length is signihcantly shorter than the model implies. Shortening this scale-length 
would cause the constraints from the rotation curve to be violated anew. We conclude 
that we can now rule out adiabatic compression of our Galaxy’s dark halo. 

Key words: dark matter - galaxies: haloes - solar neighbourhood - Galaxy: disc - 
Galaxy: fundamental parameters - Galaxy: halo 


1 INTRODUCTION 

Over the coming decade models of our Galaxy will play a 
crucial role in the extraction of physical understanding from 
the data that are currently being collected by a series of 
large surveys (e.g. Binney 2011). The models of choice are 
fully dynamical in the sense that they take full account of 
the equations of motion of the Galaxy’s constituent particles 
and of the laws of gravitation. 

An important class of such models represent the Galaxy 
by a large number of particles. The model is then specified by 
the phase-space coordinates, masses and other properties of 
the particles. Specifications of this kind are highly redundant 
in the sense that on integrating the equations of motion for a 
dynamically insignificant time, such as 10 Myr, all the phase- 
space coordinates change while the model remains the same. 
This redundancy is associated with great complexity in the 
sense that millions of particles are required to achieve even 
mediocre spatial resolution in the model, so (a) at least 10 
million real numbers are required to specify a model, and 
(b) given these numbers it is a totally non-trivial exercise 
to characterise the model in a meaningful way. Moreover, as 
a consequence of this complexity, the models are typically 
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computationally expensive to produce and it is unlikely to 
prove feasible to fit such a model to the wide range of very 
detailed data for our Galaxy that are becoming available. 

Models of a very different class are specified by a dis¬ 
tribution function (df) /a(J) for each population a of con¬ 
stituent particles. Each df is an analytic function of three 
constants of orbital motion, which it is convenient to choose 
to be the actions Ji, and specifies the probability density of 
particles of type a in the three-dimensional region of phase 
space associated with a given set J of action values. The le¬ 
gitimacy of assuming that dfs depend only on J is assured 
by Jeans’ theorem (Binney & Tremaine 2008). The model 
is specified by giving the values of the parameters p that 
appear fa alongside the actions. Hence the model is speci¬ 
fied by the values of < 100 numbers, rather than millions, 
and it is computationally feasible to adjust the parameters 
to optimise the agreement between the model and the data. 
Moreover, the parameters can be chosen to be approximately 
equal to intuitive properties of the population, such as its 
radial scale length, radial and vertical velocity dispersions 
in the solar neighbourhood, etc. 

In a series of papers Binney (2010, 2012); Binney et al. 
(2014); PifH et al. (2014) we have demonstrated the ef¬ 
fectiveness of a particular form of df for disc stars. Bin¬ 
ney (2012) fitted this df to the kinematics of stars in the 
Geneva-Gopenhagen Survey (Nordstrom et al. 2004; Holm- 
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berg et al. 2007; Casagrande et al. 2011, hereafter GCS), 
which lie within 120 pc of the Sun. The resulting model 
was then used to predict the kinematics of stars in the RA- 
dial Velocity Experiment (Steinmetz et al. 2006; Kordopatis 
et al. 2013, hereafter RAVE), which extend to ~ 2.5 kpc of 
the Sun (Binney et al. 2014). The ability of the model to 
reproduce data very different from that used in the model’s 
specification is remarkable in view of the limited extent to 
which the thick disc contributes to the GCS data. PifH et al. 
(2014, hereafter P14) showed that minor adjustments to the 
parameter values in Binney (2012) enabled the model to 
fit simultaneously the RAVE data and data from the Sloan 
Digital Sky Survey in Juric et al. (2008, hereafter JOS), and 
yielded the currently strongest constraints on the mass of 
dark matter within the solar radius, Ro- 

Fundamentally stars form a multi-dimensional contin¬ 
uum spanned by mass, age, and chemistry. Models need to 
take cognisance of the masses, ages and chemical composi¬ 
tions of stars because these intrinsic variables determine a 
star’s absolute magnitudes in the various observing bands, 
and thus the probability that a star will be included in any 
given data set. The work summarised above handled these 
complexities very crudely by supposing that our Galaxy 
comprises a continuum of stellar populations that differ only 
in age. Even within this rough approximation the analysis 
was unnecessarily crude. Sanders & Binney (2015) moved 
the use of analytic dfs to a new level of sophistication by 
(a) introducing metallicity [Fe/H] as an argument of the 
DF of disc stars, and (b) using isochrones to compute the 
number of stars the model predicts will be observed in given 
surveys. Sanders & Binney (2015) fitted their “extended df” 
for the disc to the GCS and compared the resulting predic¬ 
tions for the more distant G dwarfs in the Sloan Extension 
for Galactic Understanding and Exploration (Yanny et al. 
2009, hereafter SEGUE). 

Although much remains to be done, as a result of 
the work just summarised, we now have a reasonable un¬ 
derstanding of the phase-space distribution of disc stars - 
the extensions that seem most urgent are to compare the 
model’s predictions with data for stars that never come close 
to the Sun, and to model the way non-stationary features 
such as the Galactic bar and spiral structure modify the 
model’s static df. 

Our Galaxy is by no means comprised of just stars: 
overall more than 90 percent of its mass is thought to be 
contributed by dark-matter particles, and from the RAVE 
data using an analytic df for the disc, P14 concluded that 
even interior to the Sun, where baryons make their largest 
contribution to the overall mass, more than half the matter 
is dark. 

In the work with analytic dfs so far mentioned, dark 
matter was included merely by adding to the gravitational 
field generated by stars and gas a hypothetical contribution 
from the dark halo. The density distribution of the dark halo 
that was used by P14 and similar studies was inspired by 
simulations of cosmological clustering that did not include 
baryons, and hence side-stepped the complex physics that 
baryons bring into play. In these simulations the density 
profiles of dark halos are well approximated by a simple 
functional form, the so-called NFW profile, regardless of the 
mass or assembly histories of individual halos (Navarro et al. 
1997). Given that we have yet to actually detect any dark- 


matter particles, and dark matter strongly dominates the 
Galaxy’s overall mass budget, real dark halos are widely 
modelled as spherical or mildly spheroidal version of the 
profile. 

The intellectual basis for this procedure is weak, how¬ 
ever, because even if a given dark halo conformed to the 
NFW profile before the bulk of the baryons flowed in and 
formed stars, it would be a remarkable coincidence if it con¬ 
formed to this profile now. Indeed, if the baryons accumu¬ 
lated over many dynamical times rather than in a few mas¬ 
sive infall events then the actions of the dark-matter parti¬ 
cles will be invariant as the Galaxy grows, and it will be the 
DF of the dark matter, not its spatial density profile, which 
will be invariant. The existence of an old thin disc favours 
the hypothesis that baryons accumulate quiescently. 

This line of argument motivated (Piffl et al. 2015, here¬ 
after PPB) to open a new chapter in the use of analytic dfs 
by specifying the dark halo through its df rather than its 
spatial density distribution. Hence the Galaxy model they 
fitted to RAVE data was based on analytic dfs for both 
the stars and the dark matter, and the spatial distributions 
of these components were determined by iteratively solv¬ 
ing for the gravitational potential that these components 
jointly generate in concert with smaller, specified contribu¬ 
tions from gas and an axisymmetrised version of the Galactic 
bulge/bar. PPB took as their point of departure the model 
P14 had fitted to RAVE and JOS data. This model pro¬ 
vided the complete df of the disc and a particular NFW 
profile. PPB adopted the disc df, and for the df of the dark 
halo used one that self-consistently generates the given NFW 
profile in the absence of the disc. PPB found that when this 
dark halo cohabits with the disc, the disc’s gravitational field 
contracts the central portion of the dark halo sufficiently to 
cause the final model to be inconsistent with the measured 
circular speed Vc at R<6kpc. 

This conflict between data and the first model of our 
Galaxy to specify the dark halo through its df, can be ad¬ 
dressed in two extreme ways: one can modify the df of either 
the disc or the dark halo, leaving the other alone. Modifica¬ 
tion of the halo’s df would imply that the infall of baryons 
was sufficiently unsteady to violate adiabatic invariance of 
the actions of halo particles. Here we ask whether this con¬ 
clusion can be avoided by leaving the functional form of the 
dark halo’s df alone and seeking a disc df that is consistent 
with the observational data when the disc cohabits with a 
dark halo of this form, and the previously assumed quanti¬ 
ties of gas and bulge stars. We will show that with current 
data it is possible to model the RAVE and SEGUE data 
successfully within constraints set by measurements of Vc 
by increasing the scale length of the disc from ~ 2.5 kpc to 
3.5 kpc, but find that the resulting model then predicts 
too little microlensing of bulge stars. 

A key achievement of this paper is the development of a 
practical technique for searching a multi-dimensional model 
space for a model that satisfies observational constraints on 
the space density and kinematics of stars in the extended so¬ 
lar neighbourhood, and constraints from gas, masers and Sgr 
A* on the Galaxy’s rotation curve. In forthcoming work we 
plan to use this technique to obtain fully dynamical models 
of our Galaxy that are consistent with all available observa¬ 
tions. 

In Section 2 we list our observational inputs. In Sec- 
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tion 3 we specify the functional forms that define our mod¬ 
els. In Section 4 we explain how we fit the data. Section 5 
describes the best-fit model. Section 6 discusses the impli¬ 
cations of the model and in Section 7 we sum up and draw 
conclusions. 


2 OBSERVATIONAL INPUTS 

We use the same observations as in P14, apart from drop¬ 
ping the constraints for the dark halo. We assume a dis¬ 
tance of the Sun to the Galactic centre (GC), 7?o, to be 8.3 
kpc (e.g. Gillessen et al. 2009; McMillan 2011; Schonrich 
2012), the distance of the Sun above the Galactic plane, 
zo, to be 14 pc (Binney et al. 1997) and the solar motion 
with respect to the local standard-of-rest (LSR), vq, to be 
(11.1,12.24, 7.25) kms“^ (Schonrich et al. 2010). 

Our most important inputs are (i) the vertical profile 
of stellar density above the Sun determined by JOS, and (ii) 
the kinematics of 180 000 giant stars in RAVE. 

2.1 Vertical density profile from SDSS 

We assume that the population from which the RAVE sam¬ 
ple is drawn is identical to that studied by JOS, who mea¬ 
sured its vertical density profile by means of a main-sequence 
colour-magnitude relation. We use the data points shown in 
the middle panel of their Figure 15, which shows results 
from M dwarf stars in the colour range 0.70 < r — i < 0.80. 
Similar to RAVE, this sample should carry only weak biases 
in metallicity and age. Rather than correcting the data for 
the effects of Malmquist bias and binarity as JOS did, P14 
imposed these effects on the model. Since this step had a 
negligible effect on their results, we omitted it. 

We decomposed the JOS density profile into contribu¬ 
tions from the disc and stellar halo. This decomposition im¬ 
plied that at {R, z) = (i?o, 0.5 kpc), the density of the stellar 
halo is 0.0056 times the density of the disc. The df of the 
stellar halo is subsequently always normalised to produce 
this ratio of densities at (7?o, 0.5 kpc). 

2.2 Kinematics from RAVE 

For the kinematics we use the stellar parameters and dis¬ 
tance estimates in the fourth RAVE data release (Kor- 
dopatis et al. 2013). We define eight spatial bins in the 
{R, z) plane. Four bins for stars inside the solar cylinder 
with i?o — 1 kpc < R < Ro and \z\ in [0,0.3],[0.3,0.6],[0.6,1.0] 
or [1,1.5] kpc. The other four bins cover the same a ranges 
but cover the regions 1 kpc outside the solar cylinder, i.e. 
Ro < R < Ro + 1 kpc. After sorting the stars into these 
bins, we compute the velocity distributions predicted by the 
DF at the mean {R, z) positions (barycentre) of the stars in 
each bin. For each bin we have a histogram for each compo¬ 
nent of velocity, so we accumulate from 24 histograms. 
Throughout this work we compute velocities in the coordi¬ 
nate system that Binney et al. (2014) found to be closely 
aligned with the velocity ellipsoid throughout the extended 
solar neighbourhood - this system is quite closely aligned 
with spherical coordinates. We denote the velocity compo¬ 
nent along the long axis of the velocity ellipsoid - pointing 


more or less towards the Galactic centre - with Vi, the az¬ 
imuthal component with VI/,, and the remaining component 
with Vj. 

The resulting model distributions cannot be directly 
compared to the observed distributions, because the latter 
are widened by errors in the velocity and parallax estimates. 
We fold the model distributions with the average velocity 
uncertainties of the bin’s stars to obtain Vbary(Vi). The dis¬ 
tortions arising from the parallax error are less straight for¬ 
ward to introduce: following Binney et al. (2014) we create 
a Monte Carlo realisation of a given df by randomly assign¬ 
ing to each star in our RAVE sample a new “true” distance 
according to its (sometimes multi-modal) distance pdf, and 
a new “true” velocity according to the model velocity dis¬ 
tribution at this position. With these new phase space co¬ 
ordinates we compute new observed line-of-sight velocities 
and proper motions. These are finally equipped with random 
observational errors. Using the original catalogue distances, 
we then compute new realistically distorted velocity distri¬ 
butions, VMc(Vi), based on the df that can be compared 
directly to the original RAVE distributions in a number of 
spatial bins. We minimise the Poisson noise in AMc(Vi) by 
choosing 100 new velocities for each star. This procedure 
is computationally expensive and the distortions vary only 
weakly for reasonable choices of the df parameters. To speed 
up the process, we store the ratio Vbary(Vi)/VMc(Vi) for a 
DF that is already a good match of the RAVE data. Exam¬ 
ples of these ratios, which are near unity in the core of the 
distribution but fall to < 0.2 in the wings, are shown in the 
lower panels of Fig. 4 of P14. These ratios are then used to 
correct all df predictions before they are compared with the 
data. 

Our model selection involves computing the correspond¬ 
ing velocity histograms at the barycentre of each bin, and 
optimising the fit between the data and these histograms 
after the latter have been modified to allow for the Impact 
of errors In the measurements of velocity and distance. 

2.3 Gas terminal velocities 

The distribution of HI and CO emission in the longitude- 
velocity plane yield a characteristic maximum (“terminal”) 
velocity for each line of sight (e.g. Binney & Merrifield 1998, 
§9.1.1). The terminal velocities are related to the circular 
speed Vc{R) by 

Uterm(0 = Uc(R) - Uc (Ro) siu Z 

= Vc{Ro sinZ) — Uc(Ro)sinZ. 

We use the terminal velocities Uterm(Z) from Malhotra 
(1995). Following Dehnen & Binney (1998) and McMillan 
(2011) we neglect data at sinZ < 0.5 in order not to be in¬ 
fluenced by the Galactic bar, and we assume that the ISM 
has a Gaussian velocity distribution of dispersion 7 kms“^. 

2.4 Maser observations 

Reid et al. (2014) presented a compilation of 103 maser ob¬ 
servations that provide precise 6D phase space information. 
Since masers are associated with young stars their motions 
should be very close to circular around the GC. We again as¬ 
sume an intrinsic velocity dispersion of 7kms~^ and no lag 
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against the circular speed (van der Kruit & Shostak 1984; 
McMillan & Binney 2010). For the likelihood computation 
we neglected 15 sources that were flagged as outliers by Reid 
et al. (2014) and also all sources at i? < 4kpc. The latter 
is again to prevent a bias by the Galactic bar. To assess 
the likelihood of a maser observation, we predict the ob¬ 
served velocities (line-of-sight velocity, proper motions) as 
functions of heliocentric distance and then integrate the re¬ 
sulting probability density along the line-of-sight. 


2.5 Proper motion of SgrA* 

Reid & Brunthaler (2004) measured the proper motion of 
the radio source SgrA* in the GC to be 

MSgrA* = —6.379 ± 0.024 masyr”^. (2) 

This source is thought to be associated with the super- 
massive black hole that sits in the gravitational centre of the 
Milky Way with a peculiar velocity below lkms“^. Hence 
this measurement reflects the solar motion with respect to 
the GG. 


3 MODEL DEFINITIONS 
3.1 The stellar disc 

The functional form of the stellar disc’s df, which is made 
up of contributions from the thick disc and each coeval co¬ 
hort of thin-disc stars, is unchanged from PPB, even though 
Sanders & Binney (2015) have strictly speaking rendered 
that form obsolete. The df segments the disc by age: the 
oldest stars are represented by a df for the “thick disc”. The 
thick disc’s df is a “quasi-isothermal” component (Binney 
& McMillan 2011). The df of such a component is 

f{Jr, P, L.) = L,), (3) 

where and fa^ are defined to be 

oy' 2 

UAJr,L,) = [1 +tanh(L./Lo)]e-''‘'”/"- (4) 

TT (Jrp Kj 

and 

= (5) 

Here Q,{LA), «:(Lz) and v{Lz) are, respectively, the circular, 
radial and vertical epicycle frequencies of the circular orbit 
with angular momentum Lz, while 

Y.{Lz) = Eoe-^‘=/^^ (6) 

where Rc{Lz) is the radius of the circular orbit, determines 

the surface density of the disc: to a moderate approxima¬ 
tion the surface density at Galactocentric distance R can be 
obtained by using for Lz in equation (6) the angular mo¬ 
mentum Lz{R) of the circular orbit with radius R. 

In equation (4) the factor containing tanh serves to 
eliminate retrograde stars; the value of Lq controls the ra¬ 
dius within which significant numbers of retrograde stars are 
found, and should be no larger than the circular angular mo¬ 
mentum at the half-light radius of the bulge/bar. Provided 
this condition is satisfied, the results for the extended solar 
neighbourhood presented here are essentially independent of 
Lq. 


The DF of the thin disc is taken to be a superposition 
of quasi-isothermal dfs, one for the stars of each age r. The 
velocity-dispersion parameters ar and az above are func¬ 
tions ar{Lz,r) and (Jz{Lz,t) of angular momentum and age. 
They control the radial and vertical velocity dispersions of 
the stars of age r and are approximately equal to them at 
Rz- Given that the scale heights of galactic discs do not 
vary strongly with radius (van der Kruit & Searle 1981), 
these quantities must increase inwards, and we assume this 
dependence on Rz is exponential. We take the growth with 
age of the velocity dispersions of a coeval cohort of thin-disc 
stars from the work of Aumer & Binney (2009). With these 
assumptions the velocity-dispersion parameters are given by 


U'r(L2,T) — U 7.0 


/ r -f n 

V Fm + Tl 


g(^0~^c)/ Rt7,r 


(^z{Lzf t ) — OzO 


/ r -I- n 

V Fm + n 




(7) 


Here azo is the approximate vertical velocity dispersion of 
local stars at age Tm — 10 Gyr, ri sets velocity dispersion 
at birth, and /3 ~ 0.33 is an index that determines how the 
velocity dispersions grow with age. The radial scale-lengths 
on which the velocity dispersions decline are Ra,i, and a 
constant scale height would be follow if Rcr,z ~ 2Rd. 

We assume that the star-formation rate in the thin disc 
has decreased exponentially with time, with characteristic 
time scale to, so the thin-disc df is 


J. /T T r \ fo ^ ° f<Zr{Jr,Lz)faz{Jz,Lz) ^ 

hhAJr,Jz,Lz) - ig(er„/io _ 1) ’ 


where ar and az depend on Lz and r through equation (7). 
We set the normalising constant Eq that appears in equation 
(6) to be the same for both discs and use for the complete 

DF 


./"disc i^Jr, Jz , Lz') — Jth.n(^Jr, Jz , L z) ~\~ Fthk./thk )Jr , Jz , Lz) , 

(9) 

where Ethk is a parameter that controls the fraction (1 -f 
74hk)~^ of stars that belong to the thick disc. The values of 
the parameters for our hnal model are given in Table 2. 

We followed PPB in imposing a lower limit of 1 kpc on 
the value of RAJ A which the epicycle frequencies k{JA) 
and are evaluated for use in the df. 


3.2 DF of the stellar halo 

As in P14, we include the contribution of a stellar halo when 
htting the kinematics of RAVE stars, which include a small 
but non-negligible population of stars that are identifiable 
as halo stars by their low or even negative values of the 
azimuthal velocity v^. Including the df of the stellar halo 
prevents the fitting routine distorting the thick disc in an 
attempt to account for the presence of halo stars in the 
sample. 

The density of the stellar halo is generally thought to 
follow a power-law in Galactocentric radius, i.e. phaio oc r““, 
with the power-law index a ~ 3.5 (e.g. Binney & Merriheld 
1998, §10.5.2). We can model such a conhguration using the 
following form of the (un-normalised) df (Posti et al. 2015) 

/(J) = 5(J)"-® exp{-[<;(J)/3,aax]n. (10) 
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where 

g(J) = Jr+ 7i|L,| + 72 J, + ji:L (11) 

with 71 = 0.937, 72 = 0.682, jioL = 200kms“^kpc and 
gmax = 2.5 X 10® kms“^ kpc. These choices of 71 and 72 
make the stellar halo approximately spherical. The RAVE 
data alone are not well suited to constraining the stellar halo, 
so we defer this exercise to a later paper (Das & Binney in 
preparation). We include the stellar halo only to prevent 
distortion of the thick disc that is fitted to the data. Onr 
complete total stellar df is 

/(Jr, Jz,Lz) = fd isc(Jr, Jz, Tz ) “j“ Jhalo/halo (Jr , Jz, Ta) 

( 12 ) 

with Thaio chosen so phaio/pdisc = 0.0056 0.5 kpc above the 
Sun to be consistent with the JOS data as explained at the 
start of Section 2. 


3.3 The dark halo 


Posti et al. (2015) found a simple df /(J) that self- 
consistently generates a spherical model that has almost 
exactly an NEW profile. This model is essentially isotropic 
near its centre, and becomes mildly radially biased beyond 
its scale radius. PPB extended this df so they could explore 
how halos with different velocity anisotropies respond to the 
infall of baryons. They reported results for three model dark 
halos, one radially biased, one tangentially biased and the 
original nearly isotropic model of Posti et al. Here we con¬ 
sider only the radially biased case, which most closely re¬ 
sembles the halos that form in simulations that include only 
dark matter. The PPB df is defined in terms of an approx¬ 
imately homogeneous function of the actions 

hO)^jJr + ^i\M + JO, (13) 

where A and B are given by equations ( 6 ) and (7) of PPB 
with 6 = 8 to ensure radial anisotropy: 


A = 4.5 -b 3.5tanh^ 
B — 4.5 — 3.5tanh^ 


f \Lz\+Jz \ 

\ Jr -b I J/z I + Jz / 

{ \Lz\+Jz \ 

\ Jr -b \Lz\ + Jz / 


(14) 


The quantities k and above are epicycle and azimuthal 
frequencies in the self-consistent, isolated spherical model 
evaluated at the radius Rc of a circular orbit with angular 
momentum 


Jtot = Jr -b \ J^\ + Jz- (15) 


We take the argument of Rc to be Jtot to make it an approx¬ 
imate function of energy, so J?c does not become small, and 
the epicycle frequencies large, for stars on eccentric and/or 
highly inclined orbits with small \J^\. 

The DF vanishes for h > hmax = 10® kms“^ kpc, and 
for h < hmax 


/(J) 


N [1 + Jo/(Jcore + h(J))]®/® 

J3 [1 + MJ)/Jo] = '® 


(16) 


where N, Jo, Jcore and /tide are constants. The normalising 
constant N determines the virial mass of the dark halo, Jo 
sets the radius of the transition between the inner and outer 
power-law segments of the system’s radial density profile. 


Table 1. Parameters for gas disc and bulge. With the exception 
of Sq and Rd, these parameters were fixed. 


Parameter 

value 

unit 

Gas disc 

So 

57.8 

Mq pc ^ 

Rd 

(stars) 


Zd 

0.04 

kpc 

-^hole 

4 

kpc 

M(oo) 

17.7 X 10® 

M0 

M{Ro) 

3.53 X 10® 

M0 

Bulge 

PO.b 

94.9 

Mq pc“^ 

rg.b 

0.075 

kpc 

rcut,b 

2.1 

kpc 

<?b 

0.5 


7b 

0 


db 

1.8 


M(oo) 

8.56 X 10® 

M0 


Jcore = 1.25 X 10~^kms~^kpc is a small number required 
to keep the central phase-space density finite and /tide is 
chosen to make / vanish for h — /imax = 10 ® kms~^ kpc. 
Appendix C of PPB explains the rationale for this choice of 
the dark halo’s df. 

Once we have obtained a self-consistent model of an iso¬ 
lated halo, we freeze the dependence of Rc on its argument, 
so while we relax <l?tot onto the potential that is jointly gen¬ 
erated by all the Galaxy’s components, the dark halo’s df 
stays exactly the same function of J. It is essential to freeze 
the function /(J) during the introduction of the discs and 
the bulge if one seeks to learn how the halo is distorted by 
the gravitational helds of its companions. 


3.4 The bulge/bar and gas disc 

Our modelling technique restricts us to axisymmetric mod¬ 
els, so we cannot use a sophisticated model of the bulge/bar. 
Moreover, the data we use are only sensitive to the bulge’s 
contribution to radial forces. Therefore we do not represent 
the bulge by a df /(J) but by a fixed axisymmetric mass dis¬ 
tribution. Following McMillan (2011) we use a model similar 
to that constructed by Bissantz & Gerhard (2002). 

The density distributions of the bulge is 

p{R, z) = exp[-(mro/rcut)^l, (17) 

m'(i + m)P ' 

where 

m(R,z) = a/ {R/roY -b [zlqroY- (18) 


Our model bulge has an axis ratio q = 0.5 and extends to 
Fcut = 2.1 kpc: Table 1 lists all the parameters. 

The gas disc is likewise represented by an axisymmetric 
distribution of matter that has density 


/ p N Eg 


+ Ifi ^hole \ 

Rd Zd R ) 


(19) 


A non-zero parameter Rhoie creates a central cavity in the 
disc. The values of the parameters are given in Table 1. 
The scale length Rd is set equal to twice the scale length 
of the stellar disc, and the surface density normalisation is 
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adjusted to maintain the ratio 13.5 : 35.5 between the gas 
and stellar surface densities at Ro that is given in Flynn 
et al. (2006). Rd and Eq are the only parameters that are 
varied: the other parameters are fixed at the values adopted 
by P14 and earlier investigators. 


4 FITTING ALGORITHM 

Since the scale action Jo of the dark halo’s df has little 
bearing on the contribution of the dark halo to forces on 
stars in the inner Galaxy (r < lOkpc), we fix Jo to a value, 
Jo = 6000kms“^ kpc, consistent with values obtained in 
preliminary tests of a procedure that did involve varying Jo. 
Fixing Jo effectively constrains the virial mass of the dark 
halo, so there is then no need to impose an explicit upper 
limit on the halo mass within 50kpc as P14 did. With Jo 
determined, the only adjustable parameter for the dark halo 
is its overall normalisation. 

The disc df has in principle 12 parameters: for each sub¬ 
disc a mass, the normalising velocity dispersions ct^o and CTjo, 
and the radial scale lengths Rd , Ra ^, Ra^ of the mass density 
and the velocity dispersions. We reduce the free parameter 
count to 9 by assuming that Rd is the same for both discs, 
and for the thin disc setting Ra^ = Ra^ = 2i?d. Since the 
radial range covered by RAVE increases with \z\, the data 
do contain information about the thick disc’s values of Ra^ 
and Rcr ^, so we let these parameters be chosen by the data. 
Hence the disc df requires 9 parameters (two masses, four 
pseudo velocity dispersions, the radial mass scale length, two 
radial scale lengths for the thick disc’s velocity dispersions). 

When fitting to the data we have in all 10 free param¬ 
eters: 9 for the stellar discs, and one for the dark halo. 

The construction of a fully self-consistent model as de¬ 
scribed by PPB is computationally quite expensive because 
for each trial potential the density must be computed at ev¬ 
ery point of an extensive grid by integrating over all three 
components of velocity. Consequently, it proved impractical 
to go through the whole procedure in each iteration of the 
fitting process with in total 10 fitting parameters. To speed 
up the process, we exploit the fact that we already know in 
some detail the final mass distribution of the stellar com¬ 
ponent, because (i) our disc df always results in a density 
distribution very close to a double-exponential and (ii) the 
star-count data from JOS set very tight constraints on the 
vertical structure of the disc. Hence for large parts of the 
fitting process we do not deal with the disc df directly, but 
use a dummy potential of a double-exponential distribution 
with the right vertical density structure. Hence our proce¬ 
dure is a follows: 

1. We choose a normalisation of the total mass of the dark 
halo. 

2. We adopt a double-exponential disc that has the right 
vertical structure (satisfies the JOS data). 

3. We find the self-consistent equilibrium of our chosen 
dark halo in the presence of the bulge, the gas disc and our 
adopted double-exponential stellar disc. 

4. By adjusting masses and radial scale(s) of the double¬ 
exponential disc (and consequently of the gas disc, which 
is tied to the stellar disc), we seek satisfaction of the con¬ 
straints on Vc{R) listed in Section 2 - these comprise the 
terminal velocities and data for the masers and Sgr A*. 


Each time the the disc parameters are changed, we relax 
the halo again in the presence of the updated disc. Thus af¬ 
ter each unsuccessful comparison with data of the rotation 
curve generated by the dark-halo df in the presence of the 
current double-exponential disc, we return to step 2. 

5. We choose a df for the stellar disc that has the scale 
lengths found in Step 4 and we adopt plausible values for 
the velocity-dispersion parameters. 

6 . Then we determine a new overall potential as follows, 
(i) We update the contribution of the disc to the density 
that generates the potential to the density obtained by inte¬ 
grating the disc df over velocities using the current estimate 
of the potential, (ii) We solve for the potential generated by 
the updated disc density and the current estimate of the 
dark-halo density, (iii) We update the contribution of the 
dark halo to the density that generates the potential to the 
density obtained by integrating the halo df over velocities 
in the current potential, (iv) We again solve for the over¬ 
all potential. Then we return to step (i) until the updates 
become insignificant. 

7. With the potential frozen at its current value, we adjust 
the velocity-dispersion parameters in the disc df to obtain 
a good fit of the model’s kinematics to the kinematics of the 
RAVE giants in eight spatial bins. At this stage we include 
the contribution of the stellar halo with its df normalised 
as described in Section 3.2. Then we return to Step 6 and 
follow it with Step 7 until updates become negligible. 

8 . We compute the residuals between the vertical stellar 
density profile of JOS and that implied by the dfs of the disc 
and stellar halo. 

9. If these residuals are unsatisfactory, we choose a new 
mass for the dark halo and return to Step 2. 


This algorithm derives from two physical principles. 
First the predicted rotation curve is sensitive to the scale 
lengths of the disc components and insensitive to the ve¬ 
locity dispersion parameters. Second the velocity-dispersion 
parameters control the kinematics of stars within each spa¬ 
tial bin with little sensitivity to the potential used. Conse¬ 
quently, even when the potential employed is far from the 
truth, the RAVE data yield values for the velocity-dispersion 
parameters that are close to the final best-fit values. Given 
the value of CTzO, the correct potential can be identified by 
comparing the density profile predicted by the df to the 
star-count data. 

Replacing the double-exponential discs by the density 
distribution generated by the disc df in step 6 changes the 
potential very little. We do not include the stellar halo in 
the sources of the gravitational potential because its mass 
should be negligible, even though the mass implied by our df 
can be significant depending on its ill-constrained behaviour 
at small actions that do not contribute to the RAVE data. 
Because the density distribution generated by the disc df is 
not exactly the same as that generating the dummy poten¬ 
tial, it is not practical to set the df normalisation such that 
we obtain the same total mass as generates the dummy po¬ 
tential. Instead we normalise such that the disc df yields the 
same radial force at the Sun as does the dummy potential. 
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Figure 1. The axis ratio q of the final dark halo as a function of 
semi-major axis length a. 

5 RESULTS 

Here we describe our best-fit model. Table 3 gives the values 
of the parameters in the final dark halo’s df, and below the 
line some derived properties. Its local density is 

Pdm(Ro,0) = 0.01307 Mo pc“® = O.fiOGeVcm"®. (20) 

In the absence of the discs and bulge, the df of the dark 
halo generates a spherical mass distribution. Fig. 1 shows the 
final dark halo’s axis ratio g as a function of semi-major axis 
length a. At a < Rq the axis ratio is fairly constant at just 
below 0.8, and then, as the influence of the disc wanes, it 
increases towards q — 1, implying spherical symmetry. The 
axis ratio q ~ 0.8 in the inner Galaxy is consistent with 
the findings of PPB, but their halo became rounder sooner 
because their disc was more compact. 

For the halo scale action Jo we find a value around 
6000 km s“^ kpc. This parameter is, however, only weakly 
constrained by our data and strongly correlated with the 
derived halo virial mass because the latter is largely de¬ 
termined by the density of the dark halo well outside Rq, 
while our data constrain the halo density at r < Ro- In fact, 
the virial mass of the dark halo can be increased without 
much impact on our data by increasing Jo so the density 
within Ro changes little. Taken in isolation, i.e. used to re¬ 
cover the dark halo’s structure before the baryons fell in, our 
dark-halo df yields an NFW profile that has scale length 
Vs — 16 kpc, virial radius R 200 ~ 223 kpc, and virial mass 

Af2oo = 1.4 X 10^^ Mq. 

The stellar discs have a scale length i?d = 3.66 kpc and 
total mass 3.6 x 10^° Mq. At Ro the stellar surface density 
is 46.3 Mq pc“^. The total baryonic mass (i.e., including the 
Bulge and the gas disc) is 6.2 x 10^° Mq. Hence we have 
a (visible) baryon fraction of 4.2 per cent. The remaining 
parameters of the disc df are given in Table 2. 

In the upper panel Fig. 2 we show the final model’s cir¬ 
cular speed Vc(R)- For comparison we added in blue Vc(R) 
for the model of P14 with q = 0.8. The curves are very 
similar except for the region R ^3 kpc in which circular or¬ 
bits are impossible on account of the Galactic bar. Gonse- 
quently, the model is unconstrained in this region. The grey 
curves in this panel show circular speeds that are gener¬ 
ated by the dark matter (full curve) and baryons (dashed 
curve). The lower panel of Fig. 2 shows the ratio of the ra- 


Table 2. Parameters of the stellar disc DF. All parameters were 
adjusted when fitting the data. 


Parameter 

value 

unit 

So 

no 

M0PC-2 

Rd 

3.66 

kpc 

Thin disc 


35.40 

km s ^ 

^z,thn 

26.00 

km s“^ 

thn 

2Rd 


^(Tz ,thn 

2Rd 


Thick disc 

0'H,thk 

52.78 

km s“^ 

(^z.thk 

53.33 

km s“^ 

R ,thk 

11.6 

kpc 

R( 7 z ,thk 

5.01 

kpc 

Fthk 

0.416 



Table 3. Parameters of the dark halo’s DF. Of these parameters, 
only JV was adjusted when fitting the data. R 200 and M 200 are 
derived parameters of the dark halo that is generated by the DF 
in the absence of other components. 


Parameter 

value 

unit 

N 

8.825 x IOI 2 

Mq 

hmax 

10® 

km s“^ kpc 

Jo 

6000 

km s“^ kpc 

'^core 

1.25 x 10-2 

km s“^ kpc 

R200 

223 

kpc 

M200 

1.4 x 10^2 

Mq 


dial forces coming from the baryons and dark matter. On 
account of adiabatic contraction, our model contains much 
more dark matter at small radii than that of P14. At the 
Sun the ratio, 0.61, is even lower than that, 0.85, of P14, 
so only 38 per cent of the radial force on the Sun is due 
to baryons. Moreover, whereas in the P14 model the con¬ 
tribution to Fu from baryons rises steeply inwards, becom¬ 
ing dominant at i? < 5 kpc, in the present model it remains 
low until i? ~ 3 kpc and the baryons are dominant only at 
J? < 2 kpc. 

Fig. 3 compares the vertical stellar density profile above 
the Sun with the star count data from J08. The model star 
density comprises that predicted by the disc df plus that 
predicted by the df of the stellar halo when normalised as 
described in Section 3.2. The agreement between model and 
data is excellent. 

Fig. 4 shows our fit to the kinematics of RAVE giants in 
the four spatial bins that lie inside Rq: the coordinates of the 
barycentre of each bin are given at the bottom of the centre 
panel of each row. The fits are excellent apart from a deficit 
of almost non-rotating stars that becomes more marked as 
one moves away from the plane. This deficit points to weak¬ 
ness in our choice of df for the stellar halo and/or the thick 
disc. 

The black curve in Fig. 5 shows the vertical gravita¬ 
tional force Kz in the solar cylinder as a function of distance 
from the Galactic plane, z. The brown and blue curves show 
force laws computed from Bovy & Tremaine (2012) and P14, 
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Figure 3. Vertical stellar density profiles in the solar annulus of our composite model (solid line), the thin disc (dashed line) and the 
thick disc (dotted line), and the stellar halo (dotted line). The red and blue error bars show the data from Juric et al. (2008). Only the 
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Figure 4. Model (green curves) and observed (red error bars) velocity distributions for RAVE giants with Rq — 1 kpc < R< Rq. Each 
row represents a slice in distance from the Galactic plane. The (R, z) coordinates of the barycentre of the stars in each slice are given 
in the middle panels. These are also the locations where the DF was evaluated. A correction described in PPB was applied to the model 
prediction to incorporate the influence of errors in the measurement of both velocities and distances. 
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Figure 2. Upper panel: Vc{R) for our final Galaxy model (black 
lines) and, for comparison, for the P14 model (blue lines). The 
solid lines show the total v^, while the dashed and dotted lines 
show the contributions of the dark halo and the baryons, respec¬ 
tively. Lower panel: the ratio of the radial forces stemming from 
the baryonic components and the dark matter halo for our model 
and the P14 model. 

respectively, the blue curve being that for the P14’s model 
with q = 0.8. Error bars show results from a number of other 
studies, including the seminal work of Kuijken & Gilmore 
(1991). Only Korchagin et al. (2003), Bienayme et al. (2006) 
and P14 constrain Kz significantly at \z\ < Ikpc. In general 
there is good agreement between the studies. 

Even though our disc is slightly less massive than that of 
P14 (3.6 rather than 3.7x 10^° Mq), it generates larger Kz at 
all a because its longer scale length implies a higher surface 
density at the Sun (46.3 rather than 37.1 Mq pc“^). In fact, 
its high local surface density ensures that at all values of 2 : 
its Kz exceeds that from other studies. 


6 DISCUSSION 

By replacing the fixed NEW dark halo used by P14 with 
a dynamically active dark halo, one increases the ratio of 
dark-halo densities, pDM(3kpc)/pDM(7?o)- Since the value 
of pdm(7?o) is strongly constrained by solar neighbourhood 
data - which leave no doubt that the disc can contribute 
only part of Vc{Ro) ~ the tendency encountered by PPB for 
a more centrally concentrated dark halo to drive Uc(3kpc) 
above the observational limits must be addressed by increas¬ 
ing the scale length of the disc. If done at fixed disc mass, 
this operation moves disc material outwards, thus counter¬ 
acting the increased central concentration of the dark halo 
on its becoming active. However, as the lower panel of Fig. 6 
shows, increasing Rd at fixed disc mass reduces the disc’s 
contribution to Vc{Ro)- Hence there is little scope for reduc¬ 
ing either Mdisc or pdm(Ro) below the values found by P14. 
The upper panel of Fig. 6 shows that as Rd is increased at 
fixed Mdisc, the local stellar surface density E(f?o) increases. 



Figure 5. Vertical force fKzl in the solar cylinder as a function 
of distance 2 from the Galactic plane. See the main text for a 
further discussion 



Figure 6. Lower panel: contribution to Uc(i?o) from a thin expo¬ 
nential disc with mass Mdisc = 3-6 x Mq and varying scale 
lengths. Upper panel, the surface density of such a disc at TJq. 

The scope for such an increase is limited by the RAVE and 
JOS data. Hence the present optimisation process leads to a 
model with a longer disc scale length but similar mass and 
Pdm(Ro)- Pressure from the constraints on t;c(3kpc) drives 
the disc to higher values of E(i?o) and Kz than are ideal 
from the perspective of the RAVE and JOS data. 

6.1 Constraints from microlensing 

Above we found a model with an adiabatically compressed 
dark halo that Is consistent with all the observational con¬ 
straints of Section 2. However, we now show that this model 
Is not consistent with data not so far considered, namely the 
measured optical depth to microlensing of bulge stars. Bin- 
ney & Evans (2001) took the local density of dark matter 
to be Pdm(Ro) = (0.0136 ± 0.007) Mq pc“® in close agree- 
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Figure 7. Mass as a function of scale length for a thin exponential 
disc that contributes a given amount, IOMqpc"^, to the local 
surface density. 

ment with our value. For any adopted value of the local 
dark-matter density, the local surface density of the bary- 
onic disc followed from the measured total surface density 
Eo = 71 ± GMqpc”^ within 1.1 kpc of the plane (Kui- 
jken & Gilmore 1991). Binney & Evans assumed that the 
baryons were arranged in such a way, as regards the el- 
lipticity of the bar and the scale height of the disc, that 
the optical depth to microlensing along the line of sight 
{I = 0, |6| ~ 4°) was maximised for any given contribu¬ 
tion to Vc- Taking the optical depth on this line of sight to 
be at least 2 x 10“® (Popowski et al. 2005), and extrapo¬ 
lating the adopted local dark-matter density inwards as a 
power law of slope —a, they showed that the inferred values 
of Vc exceeded those implied by measurements of tangent 
velocities unless a was significantly smaller than unity or 
Pdm(7?o) < 0.007 Mqpc”^. Our adiabatically compressed 
dark halo has Pdm(7?o) = 0.014 and a > 1, so it comfort¬ 
ably violates these constraints. 

6.2 Are long disc scale lengths viable? 

To simplify the fitting process, we set the scale lengths of 
the thin and thick discs equal, and obtained Rd = 3.66 kpc. 
This value is in good agreement with that obtained for the 
thin disc from SEGUE data by Bovy et al. (2012). It is also 
consistent with the long scale length (Rd = 3.6 kpc) JOS 
measured for the thick disc. 

However, the thick disc is now generally distinguished 
from the thin disc by high values of the abundance of a ele¬ 
ments such as Mg relative to Fe at a given value of [Fe/H]. 
When the thick disc is defined thus, current data suggest 
that the thick disc has a shorter scale length than the thin 
disc (Hayden et al. 2015; Recio-Blanco et al. 2014; Nissen & 
Schuster 2014), contrary to the finding of JOS. If the thick 
disc does have a short scale length, this will have a sig¬ 
nificant impact on the viability of the current model, be¬ 
cause the model assigns significant mass to the thick disc 
(1.06 X 10^° Mq). The vertical density profile of JOS essen¬ 
tially fixes the contribution of the thick disc to the local sur¬ 
face density. Fig. 7 shows as a function of scale length the 
mass of an exponential disc that contributes 10 Mq pc“^ to 
the local surface density. It shows that shortening the scale 


length of the thick disc will make it more massive. Its peak 
contribution to Vc will be increased by both this increase in 
its mass and the shortening of its scale length. Moreover, 
any increase in the central density of the thick disc will fur¬ 
ther compress the dark halo and further increase Vc- Hence, 
it seems unlikely that there is significant scope for reduc¬ 
ing the scale length of the thick disc without violating the 
constraints of Section 2. 

There are indications that the chemically defined, low- 
a abundance, thin disc is flared (Minchev et al. 2015). This 
being so, our geometrically thick disc with a long scale 
length may be consistent with a short scale length for the 
a-enhanced thick disc. Gonsequently, the requirement for a 
large scale length is not as serious an issue as the require¬ 
ment for a higher optical depth to microlensing. 


6.3 Selection functions 

An aspect of the present work which should be improved, is 
the neglect of selection functions. Any survey captures only 
a fraction of the stars that exist within any volume. The 
fraction that is captured varies with distance, and with in¬ 
trinsic stellar properties such as mass, chemistry and age, 
that determine the star’s absolute magnitudes in the rele¬ 
vant wavebands. The stellar densities from JOS on which we 
have relied take fully into account the selection function of 
the SDSS photometric survey. In comparing the kinematics 
predicted by models with the RAVE data we should prop¬ 
erly have considered biases towards younger or older stars, 
since the age of a star affects its likely kinematics. Sanders 
& Binney (2015) explored the impact of accounting for age 
biases on the best-fitting parameters of dfs of the present 
type in the context of GCS data. This exercise has not yet 
been performed for RAVE data, but it should be. We doubt 
that proper treatment of age biases would materially affect 
the conclusions we have reached here. 


7 CONCLUSIONS 

We have extended the approach of PPB from the construc¬ 
tion of a single Galaxy model in which both the stellar disc 
and the dark halo are represented by distribution functions 
in the the self-consistently generated gravitational poten¬ 
tial to a systematic search through a multi-dimensional pa¬ 
rameter space of such models. As in the model of PPB the 
dark halo has the structure expected if the baryons fell in 
smoothly, so the dark matter was adiabatically compressed. 
Gonsequently, the dark halo is more centrally concentrated 
than the NFW profile, which arises in simulations of cosmo¬ 
logical clustering of dark matter only. 

Whereas the model of PPB violated the observational 
constraints on Vc{R) a,t R < Rq, our search of parameter 
space has identified a model that is consistent with these 
constraints. The secret to achieving consistency is an in¬ 
crease in the scale length of the disc at essentially fixed stel¬ 
lar mass. The derived disc scale length, Rd = 3.66 kpc, may 
be acceptable for the thin disc, but it is probably unaccept¬ 
ably large for the thick disc. However, our work suggests 
that a model with an adiabatically compressed dark halo in 
which the thick disc has a realistically short scale length can- 
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not be made consistent with the constraints on Vc at small 
radii. 

Since dark matter does not cause microlensing and the 
model assigns to the dark halo much of the overall density 
at small radii, it violates by some margin the constraints on 
the dark halo that Binney & Evans (2001) deduced from the 
microlensing data. 

Hence by ruling out adiabatic compression this study 
furnishes compelling evidence that there has been a signif¬ 
icant transfer of energy from the baryons to the dark halo. 
Such a transfer has been proposed by many authors because 
it can arise through drag on the bar (Tremaine & Weinberg 
1984; Sellwood 2008), drag on molecular clouds (Nipoti & 
Binney 2015), energy injection by supernovae (Mashchenko 
et al. 2008), and by cosmic infall (Johansson et al. 2009; 
Goerdt et al. 2010). 

Given that we can now exclude adiabatic compression 
of the dark halo, the natural next step is to impose priors 
on the scale lengths of the thin and thick discs derived from 
spectroscopic studies of Galactic chemistry, and to add the 
microlensing data to the constraints listed in Section 2. Then 
one should use the fitting procedure developed here to build 
self-consistent models using dfs for the dark halo that im¬ 
ply increasing extents of central heating of dark matter by 
baryons. In this way it should be possible to place a lower 
limit on extent of dark-matter heating. We hope to publish 
the results of such a study in the near future. 
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